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ABSTRACT 

We present hydrodynamical simulations of the photoevaporation of a cloud with large-scale density gradients, 
giving rise to an ionized, photoevaporation flow. The flow is found to be approximately steady during the large 
part of its evolution, during which it can resemble a "champagne flow" or a "globule flow" depending on the 
curvature of the ionization front. The distance from source to ionization front and the front curvature uniquely 
determine the structure of the flow, with the curvature depending on the steepness of the lateral density gradient 
in the neutral cloud. 

We compare these simulations with both new and existing observations of the Orion nebula and find that a 
model with a mildly convex ionization front can reproduce the profiles of emission measure, electron density, 
and mean line velocity for a variety of emitting ions on scales of 10 17 -10 18 cm. The principal failure of our 
model is that we cannot explain the large observed widths of the [O i] 6300 A line that forms at the ionization 
front. 

Subject headings: H II regions, Hydrodynamics, ISM: individual (Orion Nebula) 



1. INTRODUCTION 

H ii regions, formed by the action of ionizing photons from 
a high-mass star on the surrounding gas, often occur in highly 
inhomogeneous environments. Most Hn regions are found 
close to high density concentrations of molecular gas and, 
for those regions that are optically visible, there is a ten- 
dency for the ionized emis sion to be blu eshifted with respect 
to the molecular emission (Israel 1978). Since optically vis- 
ible Hn regions must be on the near side of the molecular 
gas, this implies that the ionized gas is flowing away from 
the molecular cloud. The archetypal and best-studied exam- 
ple of such a "blister-type" Hn region is the nearby Orion 
nebula, which is ionized by the Q7 V star, 6 l Ori C (for re - 
cent reviews, see lO'DelJ20'01atlFerlandl200UIO'DelJ2001bl) . 
This compact Hn region, which lies in front of the molec- 
ular cloud OMC 1, shows a bright core of ionized emis- 
sion, approximately 10 18 cm in diameter, surrounded by a 
fainter halo that extends o ut to distances of order 10 19 cm 
( Subrahma nvan et alJ2001l) . Surprisingly, the thickness along 
the line of sight of the layer responsible for the bulk of the ion- 
ized emission in the core of the nebula is only ~ 2 x 10 17 cm 
(e.g.. lBaldwin et alJl99ltlWen & Q'Dellll995l) . which is sig- 
nificantly smaller than both the lateral extent of the core and 
the deduced distance between the ionizing star and the front, 
which is also of order 10 18 cm. Three-dimensional recon- 



struction of the shape of the nebula's principal ionization front 
SWen & 0'Dellll 995) shows the front to have a saddle-shaped 
structure, with a slight concavity in the NS direc tion and a 
slight convexity in the EW direction (see lO ' DeHl200 1 bl Fig- 
ure 3). In both cases the radius of curvature is large com- 
pared with star-front distance, although there are also many 
smaller scale irregularities in the front, with the most promi- 
nent of these being associated with the "Bright Bar" in the SE 
region of the nebula. Extensive spectroscopic studies of the 
Orion nebula have been carried out to determine the velocity 
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structure in the ionize d gas JKaleHll967l iBalick et al.lll974t 
[ Pankonin etail Il979t lOTtell et all 119931 iHennev & O'Delll 
119991 lO'Dell et all 120011 iDoi et alJEOOl . Although these 
show a wealth of fine-scale structure, the fundamental result 
is that lines from the more highly ionized species such as Ha 
and [Om] are blue-shifted by approximately lOkms -1 with 
respect to the molecular gas of OMC-1, which lies behind 
the nebula. Lines from intermediate-ionization species such 
as [S n] and [N n] being found at intermediate velocities. All 
these considerations lead to a basic model of the nebula as 
a thin layer of ionized gas that is strea ming away from th e 
molecular cloud, as was first proposed by Zu ckermanllll973l) . 

Despite, or perhaps because of, the mass of observational 
data that has been accumulated on the Orion nebula, no at- 
tempt has been made to develop a self-consistent dynamical 
model of the ionized flow since the early work of Vand ervoortl 
Most modern models of the nebula have con- 
centrated great effort on detail s of the atomic physics an d ra- 
diative transfer iBaldwin et alJll991t iRubin et alii 19911) and, 
while these have been quite successful in explaining the ion- 
ization structure of the nebula, they implicitly assume a static 
configuration and are forced to adopt ad hoc density con- 
figurations that have no physical basis. The importance of 
photoevaporation flows has been reco gnized in the contex t 
of the model i ng of other Hii regions iHester e t alJ (I19 9"6T): 
IScowen et ail (119981) : ISankrit & Hesterl (120001) : iMoore et all 
(2002) but again no attempt has been made to treat the dy- 
namics in a consistent manner. More recently, steady-state 
dynamic models of weak-D ionization front s have been de - 
veloped and applied to the Orion nebula ( Hennev et al. 2005), 
with some success in explaining the structure of the Bright 
Bar, where the ionization front is believed to be viewed edge- 
on. However, the one-dimensional, plane-parallel nature of 
these models renders them incapable of capturing the global 
geometry of the nebula and indeed they completely fail to re- 
produce the observed correlation between line velocity and 
ionization potential. 

The evolution of an Hi i region in an inho mogeneous 
medium was first studied bv iTenorio-Taglel fl 979 ). who con- 
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sidered the "break-out" of a photoionized region from the 
edge of a dense cloud, such as would occur if a high-mass 
star formed inside, but close to the edge of, a molecular 
cloud. Once the ionization front reaches the edge of the 
cloud, it rapidly propagates through the low-density inter- 
cloud medium, followed by a strong shock driven by the 
higher pressure ionized cloud material. At the same time, 
a rarefaction wave propagates back into the ionized cloud, 
giving rise to a "champagne flow" in which the Hn region 
is ionization-bounded on the high density side and density- 
bounded on the low-density side and the ionized gas flows 
down the density gradient. This work has been extended in 
many subsequent papers using both one-dimensional and two- 
dimensional numerical simulations ( Bodenheimer eLalil979t 
Bediin & Tenor ioTaglel 1 981: Yor ke et alJl9 82: Franc o et al.l 
1989; Comeron 119971) . which generally concentrate on non- 
steady aspects of the flow. 

Another scenario that has been extensively studied is the 
ionization of a ne u tral globule by a n external radiation field 
( Pottasch 1956; Dvson 1968; Tenorio-Tagle & Bediin 
1 19811 iBediin & Tenorio-Tagld fl984t iBertoldil 119891: 
Bertoldi & McKee 1990), to form what is variously re- 
ferred to as an ionized photoevaporation or photoablation 
flow. In this case, after an initial transient phase in which the 
globule is compressed, the ionization front propagates very 
slowly into the globule and the ionized flow is steady to a 
good approximation. Although originally developed in the 
context of bright rims and globules in H n regions, similar 
models have more recently been appli ed to the photoevapora- 
tion of circumstellar accretion disks (Johnstone et al. 1998; 
iHennev & Arthu rl 119981) and com e tary knots in plan etary 
nebulae ( Lope z-Martin et alJ 1200 ll lO'Dell et all 120001) . In 
all these cases, the ionization front has positive (convex) 
curvature, and the ionized gas accelerates away from the 
ionization front in a transonic, spherically divergent flow. 

In a strictly plane-parallel geometry, a time-stationary 
champagne flow cannot exist. This is because, in the absence 
of body forces, the acceleration of a steady plane-parallel flow 
is proportional to the gradient of the sound speed, which is 
roughly constant in the photoionized gas. However, it has 
recently been shown llHennevl|2003l) that a steady, ionized, 
isothermal, transonic, divergent flow can exist even if the ion- 
ization front is flat, so long as the ionizing source is at a finite 
distance from the front. In such a case, the divergence of the 
radiation field breaks the plane-parallel geometry and induces 
a (weaker) divergence in the ionized flow, allowing it to ac- 
celerate (an expanded and corrected version of this derivation 
is presented in Appendix lAl below) . It is natural to speculate 
that a steady flow may also be possible even if the front has 
negative (concave) curvature, so long as the radius of curva- 
ture of the front is greater than its distance from the ionizing 
source. By this means, the champagne flow problem can be 
effectively reduced to a modified version of the globule flow 
problem. 

In order to investigate this possibility, we here present hy- 
drodynamical simulations of the photoionization of a neutral 
cloud with large-scale density gradients and develop a tax- 
onomy of the resulting H n region flow structure and dynam- 
ics. We concentrate on the quasi-stationary stage of the H n 
region evolution, in which the structure changes slowly com- 
pared with the dynamical timescale of the ionized flow, with 
the result that the flow is approximately steady in the frame of 
the front. Depending on the curvature of the i onization fro nt, 
our models can resemble both "blister flows' ' Jlsrae1H97 8) if 



the cu rvature is negative (concave), or "bright rims" JPottaschl 
(19561) - if the curvature is positive (convex). The algorithm, 
initial conditions, and parameters of these numerical simula- 
tions are described in §[2]and their results are set out in §|3] In 
§|4]we derive empirical parameters for the Orion nebula that 
will allow us to make meaningful comparisons with our mod- 
els. These comparisons are carried out in §|5] In §[6] we dis- 
cuss the implications of our model comparisons to Orion and 
outline some shortcomings of the present simulations that will 
be addressed in future work. In §0we draw the conclusions 
of our paper. Finally, we present some more technical mate- 
rial in two appendices. The first derives an approximate ana- 
lytic model for the flow from a plane ionization front, while 
the second investigates the question of whether the ionization 
fronts in our simulations are D-critical. 

2. NUMERICAL SIMULATIONS 

In this section we present time-dependent numerical hydro- 
dynamic simulations of the photoionization of neutral gas by 
a point source offset from a plane density interface. 

2.1. Implementation Details 

The Euler equations are solved in two-dimensional cylin- 
drical geometry, using a second-order, finite-volume, hybrid 
scheme, in which alternate timesteps are calculated by the Go- 
dunov method and the Van Leer flux- vector splitting method 
(Godunov 1959: Ivan Leerl[l982l) . 2 The radiative transfer of 
the ionizing radiati on is carried out by the method of short 
characteristics (Raga et al. 1999), with the only source of ion- 
izing opacity being photoelectric absorption by neutral hy- 
drogen. Only one frequency of radiation is included but 
radiation-hardening is accounted for by varying the effective 
photoabsorption cro ss-section as a function of optical depth 
(Henn ev et alJl2005l) . The diffuse field is treated by the on- 
the-spot approximation. Dust opacity is not explicitly in- 
cluded in our calculations but is treated to zeroth order by 
merely reducing the stellar ionizing luminosity to account for 
absorption by grains. As shown in AppendixlAl this is a rather 
good approximation at least in the case of a flow from a plane 
ionization front. 

Ionization and recombination of hydrogen is calculated by 
exact integration of the relevant equations over a numerical 
timestep. Heating and cooling of the gas are calculated by a 
second-order explicit scheme. The only heating process con- 
sidered is the photoionization of hydrogen, with the energy 
deposited per ionization increasing with optical depth accord- 
ing to the hardening of the radiation field. The most important 
cooling processes included are collisional excitation of Ly- 
man alpha, collisional excitation of neutral and ionized metal 
lines (assuming standard ISM abundances), and hydrogen re- 
combination. Other cooling processes included in the code, 
such as collisional ionization and bremsstrahlung, are only 
important at higher temperatures than are found in the current 
simulations. 

2.2. Initial Conditions 

At the beginning of the simulation, the ionizing source is 
placed at r — 0, z — Zo, and the computational domain is filled 
with neutral gas with a separable density distribution: 

p(r,z)=p F(r)G(z). (1) 

2 The velocity shear in the Van Leer steps helps avoid the "carbuncle insta- 
bility", which occurs in pure Godunov schemes when slow shocks propagate 
parallel to the grid lines i Ouirk 1994; San ders et alll99a) . 
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For comparison with the analytic model of Appendix lAl we 
wish to contrive a neutral density distribution that will ensure 
that the ionization front stays approximately fiat throughout 
its evolution, although we will also investigate the flow from 
curved fronts. 

The initial flatness of the ionization front is guaranteed by 
choosing G(z) to be a step function between a very low value 
(that will be ionized in a recombination time) and a very high 
value (that will trap the ionization front at the interface). In 
practice, we choose the function 



log 10 G(z) = 6 tanh((z - Z\)lh x ), 



(2) 



where h\ is the thickness of the interface, which occurs at a 
height z\ (> Zo). The density contrast between the two sides 
is 10 26 . 

In order for the ionization front to remain flat as it propa- 
gates into the dense neutral gas above the interface, we let the 
initial density fall off gradually in the radial direction accord- 
ing to the function 



F{r) = 



1 + 



Zl - Zo 



(3) 



where a is an adjustable parameter controlling the steepness 
of the distribution. The analytic model of Appendix [A] sug- 
gests that a — 1 should ensure that all parts of the ionization 
front propagate at the same speed. Choosing smaller values 
of a should mean that the front propagates fastest along the 
symmetry axis and so the front becomes concave with time. 
Conversely, larger values of a will cause the front to prop- 
agate faster away from the axis and so become convex with 
time. 

A characteristic temperature, Tq, is chosen, which corre- 
sponds to gas at a density of po. The temperature distribution 
throughout the grid is then set so as to give a constant pressure 
where p > po and a consta nt tempera ture where p < po? 

In this paper we follow Shu (1992) in denoting gas veloc- 
ities in the frame of reference of the ionizing star by u, gas 
velocities in the frame of reference of the ionization front by 
v, and pattern speeds of ionization, heating, and shock fronts 
by U. 

2.3. Model Parameters 

We present three different models, which illustrate the effect 
of different radial density distributions in the neutral gas, with 
a = (Model A), a = 1 (Model B), and a = 2 (Model C). 
All other parameters are identical between the three models 
and are chosen to be representative of a compact H n region 
such as the Orion nebula. The simulation is carried out on a 
512x512 grid, with sides of length 1 pc. The source has an 
ionizing luminosity g H = 10 4 V, with a 40, 000 K black- 
body spectrum, and is placed half-way up the simulation grid 
on the z-axis. The characteristic density for the simulation is 
set at po = 10 4 m crrT 3 , where m is the mean mass per nucleon 
(assumed to be 1.3mn). The separation, z\ - Zo, between the 
source and the density interface is chosen so that the analytic 
steady-state model of Appendix lAl below would give a peak 
ionized density of po/ra = 10 4 cirT 3 in the photoevaporation 
flow. This yields zi-zo = 4.8xl0 17 cm = 0.155 pc. The initial 

3 We require a constant pressure in the high-density, undisturbed gas so 
that it remains motionless throughout the simulation. The low density gas is 
designed to be immediately ionized and then swiftly swept from the grid by 
the photoevaporation flow. A constant temperature ensures that this happens 
as smoothly as possible. 



density contrast across the interface is set at 6 = 2, yielding 
n = 10 6 cm~ 3 on the high-density side and n = 10 2 cm -3 on 
the low-density side. The thickness of the interface is set at 
h\ = 0.3(zi - zo) and the characteristic temperature is set at 
r = 600K. 4 

In order to make the simulations feasible, it is necessary to 
use a hydrogen photoionization cross-section that is consid- 
erably smaller than the physically correct value, which would 
be roughly <x = 3 x 10 ~ 18 cm 2 for unharden ed radiation from 
a 40, 000 K black body jHennev et alJl2005l Appendix A). In 
these models, we instead use a value of o~o — 5 x 10~ 20 cm 2 , 
which is 60 times smaller. This is necessary in order to prop- 
erly resolve the deep part of the ionization front (at t - 10) 
where the initial heating of the photoevaporation flow oc- 
curs. Failure to resolve this region results in spurious large- 
magnitude oscillations in the flow. Using a smaller value of 
o"o increases the width of the ionization front but has a neg- 
ligible effect on the kinematic behavior of the gas once it is 
fully ionized. We have studied the effect of increasing or de- 
creasing cr and find that our results are generally not affected, 
with the exception of some specific points noted below. 

Although each model is calculated for specific values of 
the length scale and ionizing luminosity, they may be safely 
rescaled to different values so long as the ionization parame- 
ter remains the same, which will guarantee that the ionization 
fractions in the model are not affected by the rescaling. This 
is because the dimensionless form of the gas-dynamical, radi- 
ation transfer and ionization equations co ntain only a h andful 
of dimensionless parameters (for example lHennev et al.l2005l 
Appendix A). The only one of these that depends on the den- 
sity and length scales is a dimensionless form of the ionization 
parameter: t» = noCo^o, where (q is the length scale. 

One can then use the Stromgren relation, (2h x n o^O' t0 
eliminate the density 5 and find that a fixed ionization param- 
eter requires a fixed ratio of (2h/A)- AH times should also be 
scaled proportionately to (. Furthermore, it is even possible 
to scale Qu and 6 independently, subject to the understanding 
that one is also thereby scaling the effective photoionization 
cross-section as cr oc ({/Qh) 1 ^ 2 - 

3. MODEL RESULTS 
3.1. Early Evolution of the Models 

The initial evolution of all three models is the same: the 
low-density gas below the interface is completely ionized on a 
recombination timescale of ~ lOOOyr. At the interface itself, 
the ionization/heating front at first penetrates to gas of density 
~ po, which is heated to 9000 K, and achieves a thermal 
pressure 30 times higher than the pressure in the cold neutral 
gas. This drives a Shockwave up the density ramp at a speed 
M s — (Pi/Pn)^ 2 Ci, where p\ is the ionized density, c; is the 
ionized sound speed, and p n is the pre-shock neutral density. 6 
The shock can be first discerned in the simulations after about 
3000 yr, at which time its speed is u s 5kms _1 , falling to 

1 km s _! by the time it reaches the top of the density ramp 
(after about 40, 000 yr). The time-evolution of the shock front 

4 Note that this implies temperatures as low as 10 K in the neutral gas. We 
do not attempt to model the heating/cooling or molecular processes in the 
cold gas. 

5 Although ionizations exceed recombinations by a small (=! 10/t.) frac- 
tion, this fraction is constant for all models with the same r» . 

6 The shock actually goes a bit faster than this equation suggests since the 
pressure at the heating front (see below) is about twice the pressure in the 
fully ionized gas. 
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Fig. 1 . — (a) Evolution with time of the axial distance from the star to the 
shock front, z s (upper line), heating front zi, (middle line), and ionized density 
peak, z p (lower line), all for Model B. (b) Evolution of z p for the different 
models: A (solid line), B (dashed line), and C (dotted line), (c) Same as b but 
for the evolution of z,h- (d) Same as b but for the evolution of z s . 



position in the three models is almost identical during this 
early stage, as can be seen in Figure[Q/. 

At the same time as the shock propagates into the neutral 
gas, the warm gas expands down the density ramp, towards 
(and beyond) the star. Within 30, 000 yr all the gas that was 
originally below the density interface has been swept off the 
grid by this flow, which accelerates up to about 20 km s . Ini- 
tially, the ionized gas flows parallel to the z-axis but as time 
goes on the streamlines begin to diverge and the flow settles 
into a quasi-steady configuration, during which the properties 
of the ionized flow and the subsequent evolution of the ion- 
ization front vary significantly between the models. 

In addition to the position of the shock front, z s , we make 
reference below to the position of the heating front, Zh, which 
is a very thin layer at the very base of the evaporation flow, 
and also to the position of the ionized density peak, z p , formed 
by the opposing gradients in total density and ionization frac- 
tion. The time evolution of these three quantities is shown in 
Figure [2 

3.2. Properties of the Quasi-steady Flow 

After the initial transient phase, the flow settles down into a 
long quasi-steady phase in which the flow properties change 
only on timescales much longer than the dynamic time for 
flow away from the ionization front. As we will show be- 
low, the flow structure, once normalized to the distance of the 
star from the front, seems to depend only on the curvature of 
the ionization front. During this quasi-steady phase the shock 
in the neutral gas is propagating at about 1 kms -1 , as dis- 
cussed in the previous section, and the post-shock gas moves 
at about three-quarters of this speed since it does does not cool 
strongly in our simulations. At the same time, the relative ve- 
locity between the ionization/heating front and the neutral gas 
immediately outside it is very low (for a D-critical front this 
would be c n /2cf =t 0.05 km s _1 , see Appendix 151. so that the 
ionization front speed is also about three-quarters of the shock 
speed, leading to a steady thickening of the shocked neutral 



layer. However, this behavior is dependent on the cooling be- 
havior of the shocked neutral gas, which we do not treat very 
realistically. Stronger cooling would cause the ionization and 
shock fronts to propagate at the same speed. 

Details of the variation along the z axis of various physical 
quantities for Model B are shown in Figure |2 Images of the 
distributions of ionized density and pressure in the models at 
a representative point in their evolution are shown in Figure[3] 
together with the flow streamlines. 

3.2.1. Structure along the z-axis 

The shock driven into the neutral gas is traced by the pres- 
sure (dotted line in Fig.|2jj, grayscale in right-hand panels of 
Fig. |3}, which peaks immediately behind the shock. The ion- 
ization front is traced by the ionization fraction (dot-dashed 
line in Fig. [2jj) and the ionized density (solid line in Fig.|2jz) 
shows a peak just inside the front. One can also see a heating 
front, best seen in the temperature (dashed line in Fig. |3Jj), 
which lies just outside the ionization front. The temperature 
shows a slight peak at the ionization front due to the radiation 
hardening. 

The solid line in Figure |2 shows the negative of the gas 
velocity (i.e. velocity back towards the ionizing star). The 
neutral gas is initially accelerated away from the star by the 
shock, which gives the neutral shell a speed of — hi kms -1 . 
In the ionization front it is accelerated back towards the star, 
quickly reaching speeds of ~ -10kms~' and then contin- 
uing a gradual acceleration to ~ -20kms~' before leaving 
the grid. The isothermal sound speed is shown by the dashed 
line in Figure |2j> and it can be seen that the flow becomes su- 
personic shortly after passing the peak in the ionized density. 
The question of whether the fronts are D-critical is addressed 
in Appendix |B] The net rate of energy deposition into the 
gas (photoionization heating minus collisional line cooling) 
is shown by the thick gray line in Figure|2j> and always shows 
two peaks. These can be seen in more detail in Figure 0] 
which shows a log-log plot of various quantities as a function 
of the optical depth to ionizing radiation, r = f{o-)n(Ho) dz. 
The deeper of the two energy deposition peaks is the heat- 
ing front, which occurs at an ionizing optical depth of r - 8, 
where the ionization fraction is around x 0.01 and the gas 
density around n a 10 5 cm 4 . In this front the gas is rapidly 
heated up to just under 10 4 K and accelerated up to 2-3 km s _1 
but remains largely neutral (note that in Fig.0]the velocity is 
now shown in the rest frame of the heating front, v). The sec- 
ond, broader peak, at t = 2-4, coincides with the ionization 
of the gas at a roughly constant temperature and its gradual 
acceleration up to the ionized sound speed. 

3.2.2. Differences between the models 

Despite the general similarities, there are some striking dif- 
ferences between the physical structure of the three models, 
as can be seen from Figure[3] In Model A (a — 0), where the 
neutral gas has a constant density in the r-direction, both the 
shock front and the ionization front are concave from the point 
of view of the star, whereas in Model B (a — 1), where the 
neutral gas density falls off with radius aymptotically as r~ 2 , 
the fronts are approximately flat. An analytic treatment of the 
flow in this special case is given in Appendix 1X1 In Model C 
(a = 2), which has an even steeper radial density profile, the 
fronts are convex from the point of view of the star. As a 
result of these changes, one sees several clear trends going 
from Model A, to B, to C — the ionized density falls off faster 



Flows from Smooth Ionization Fronts 



5 




CD 

B 

u 



-0.5 



-0.25 





z (parsecs) 



0.25 



0.5 



20 



15 - 



10 - 



5 - 




-0.5 



-0.25 





z (parsecs) 



0.25 



0.5 



Fig. 2. — Axial profiles of physical variables in Model B at an age of 78, 000 yr. (a) Ionized hydrogen density r\\ (solid line, units of 10 cm ), gas temperature 
T (dashed line, units of 10 4 K), pressure P g!is (dotted line, units of 10~ 7 dynecirT 2 ), and ionization fraction (dot-dashed line), (b) Negative of the gas velocity in 
the reference frame of the star, u (solid line), and isothermal sound speed (dashed line, both in km s~ 1 ), plus net heating rate (thick gray line, arbitrary units). 



in the direction towards the star; the ionized flow accelerates 
more sharply; the flow streamlines show greater divergence; 
the pressure driving the shock is greater; and, the ionization 
front progresses faster at late times. 

The lateral structure (Figure [5} is roughly similar for all 3 
models, except for the electron density, which falls off more 
rapidly in Model C, and the flow thickness close to the axis, 
which is significantly larger for Model A. 

3.3. Secular Evolution of the Flow in the Quasi-Steady 
Regime 

Figure [5] shows how the effective thickness (panel a) and 
sonic thickness (panel b) of the flows vary with time for the 
different models. It can be seen that once the quasi-steady 
regime is reached (t > 5 x 10 4 yr), the flow thickness increases 
with time for Model A, but decreases with time for Model C, 
with Model B showing a gentler decrease. 

The curvature k of each surface z{r) was found by fitting 
a parabola, z — Zo + 0.5at 2 , to the surface for r < 0.25 pc. 
Thus kl -1 is equal to the radius of curvature of the surface, 
with positive k corresponding to a convex surface that curves 
away from the ionizing star and negative k corresponding to 
a concave surface that curves towards the ionizing star. Fig- 
ures and d show the evolution with time of the dimension- 
less curvature, kzo, calculated for the heating front and for the 
maximum in the ionized density, respectively. 

Figure graphs the dimensionless effective thickness 
against dimensionless curvature for all models and all times. 
It can be seen that the vast majority of points approximately 
follow a single curve for the three models. 7 This is consis- 
tent with the idea that the flow thickness is controlled by the 
curvature of the ionization front. When the front is concave 
(/Ch < 0), the flow is not very divergent until it has passed 

7 The points that lie below this curve, at small values of k, correspond to 
the initial non-steady evolution at the beginning of each model run. 



the ionizing star, leading to a flow thickness of order the star- 
front distance. When one moves to flat and convex fronts, 
the flow divergence becomes more pronounced and its thick- 
ness becomes considerably less than the star-front distance. 
In the limit of large convex curvature («hZh,o 2> 1)> then one 
wo uld expect the flow to rese mble that from a cometary glob- 
ule (Bertoldi & McKee 1990), with the flow thickness propor- 
tional to the radius of curvature: h = cj/k, where to 0.12. 
This curve is shown by the thick gray line in Figure Q and 
it can be seen that our models are still some way from this 
regime. 

3.4. Optical Line Emission 

We have calculated the emissivity for three different types 
of line, representative of the strongest optical lines emitted 
by Hn regions. The first is a recombination line, such as the 
hydrogen Balmer lines, with emissivity 



ty-ecO) « n e n j+1 T 



(4) 



where nj+i is the number density of the recombining ion, 
which we take to be proportional to n\. The second and third 
are both forbidden metal lines, excited by electron collisions, 
with emissivity 



?7col(*) ' 



1 + B co in e T 



-1/2 



,-EjkT 



(5) 



where E is the excitation energy of the upper level, which 
we fix at /ic/6500 A, and Z? co i is the collisional de-excitation 
coefficient, chosen to give a critical density of 1000 cm -3 at 
10, 000 K. The second line is emitted by a neutral metal 
with tij oc n n , whereas the third line comes from an ionized 
metal with rij oc «j. Although the emissivities that we use are 
generic, they are each "inspired" by a particular optical emis- 
sion line. In particular, the recombination line represents a 
hydrogen Balmer line such as Ha, the collisional neutral line 
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Fig. 3. — Snapshot images of the three models (top to bottom: A, B, C) at a time of 78, 000 yr. Left column shows ionized hydrogen density as a linear negative 
grayscale, where black represents a density of 6600 cirT 3 . Right column shows gas pressure (black = 6.1 X 10~ 8 dyne cirT 2 ) with superimposed streamlines 
whose darkness indicates magnitude of the gas velocity (black = 20.1 kms~'). All grayscales are common between the models. The position of the ionizing star 
is indicated by a filled circle in each panel. 



represents [O i] 6300 A, and the collisional ionized line repre- 
sents a combination of [N n] 6584 A and [O m] 5007 A (since 
we do not consider the ionization of helium in our models, we 
cannot distinguish between singly and doubly ionized metals). 

Figure shows the emissivity of each line as a function of 
Z for each model, with arbitrary normalization. It can be seen 
that the neutral collisional line is confined to a thin region 
around the ionization front. This is because its emissivity con- 
tains a factor of n e n n oc x(l - x), which peaks at x — 0.5. The 
other two lines show emissivity profiles that are broadly sim- 
ilar to one another since they are both proportional to The 
differences between these two are mainly due to their different 
temperature dependences: the recombination line emissivity 
is relatively stronger at lower temperatures near the ionization 
front, while the collisional line is relatively stronger at higher 
temperatures, farther out in the flow. This behavior is further 
accentuated by the collisional deexcitation of the metal line, 
which reduces its strength at the higher densities found close 
to the ionization front. 

We have also calculated the spectroscopic emission line 
profiles that would be seen by an observer looking along the 



z-axis. We assume that the radiative transfer can be calculated 
in the optically thin limit, in which case the emergent intensity 
profile is 



/(H) 



f 

Jo 



7](z) exp 



dz, 



(6) 



2A 2 (z) 

where u is the observed Doppler shift and A(z) is the thermal 
width at each point in the structure, which depends on the 
local sound speed and on the atomic weight, A, of the emitting 
species. The results are shown in Figure [5] assuming A = 1 
for the recombination line, A — 14 for the ionized metal line, 
and A = 16 for the neutral metal line. 

The neutral metal line shows only a small blueshift (^ 
2kms~ 1 ) of its peak with respect to the quiescent gas, but 
also has an extended blue wing. Its profile is almost indis- 
tinguishable between the different models. The ionized metal 
line shows the greatest variation between models but all are 
qualitatively similar with a large net blueshift (12-17 km s" 1 ) 
and a FWHM of lOkms , which is roughly double the 
thermal width. In cases in which the viewing angle has been 
determined by independent means, this line may be used to 
discriminate kinematically between the different models. The 



Flows from Smooth Ionization Fronts 



7 



0.01 



0.001 



— I — I — I — I I I I I 




1 10 
Optical depth, r 



100 



Fig. 4. — Structure of the ionization front in Model B at an age of 78, 000 yr. 
Gas velocity with respect to the heating front, v, in kms~' (solid line); ion- 
ization fraction, x (dashed line); isothermal sound speed, c, in kms~' (dotted 
line); temperature, r/10 4 K (dot-dashed line); net energy deposition rate in 
arbitrary units (thick gray line). 
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Fig. 5. — Lateral variation of various model parameters at an age of 
78, 000 yr. All lengths are normalized to zh.o, the distance between the star 
and the heating front along the axis, and are plotted against cylindrical radius. 
The different models are shown by: A (solid line), B (dashed line), and C 
(dotted line), (a) Peak ionized density at each radius, normalized to value on 
the axis, (b) Total ionized emission measure, integrated along the z-direction: 
EM = Jn? dz. (c) Sonic thickness, /i SO mc, defined as the z-distance between 
the heating front and the sonic point, (d) Effective recombination thickness, 
h e g = EM/Hp. (e) Position of the ionized density peak, z p . (f) Position of the 
heating front, z h . 
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Fig. 6. — Evolution with time of the flow thickness along the z-axis and cur- 
vature, normalised to the distance, zh.o. between the star and the heating front. 
(a) Effective recombination thickness, (b) Sonic thickness, (c) Curvature of 
heating front, (d) Curvature of ionized density ridge. The different models 
are shown by: A (solid line), B (dashed line), and C (dotted line). The thin 
horizontal line in the bottom two panels indicates zero curvature and is the 
dividing line between convex and concave surfaces. 
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versus heating front curvature. Lines represent models A (solid), B (dashed), 
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Fig. 8. — Axial profiles of emissivity in different emission lines at an age 
of 78, 000 yr. The different models are shown by: A (solid line), B (dashed 
line), and C (dotted line), (a) Collisionally excited neutral metal line, (b) Col- 
lisionally excited ionized metal line, (c) Recombination line. 
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Fig. 9. — Simulated emission line profiles for a face-on viewing angle at 
an age of 78, 000 yr. The different models are shown by: A (solid line), B 
(dashed line), and C (dotted line), (a) Collisionally excited neutral metal line 
with atomic weight A = 16. (b) Collisionally excited ionized metal line with 
A = 14. (c) Recombination line with A = 1. 



recombination line is also blueshifted (by about lOkms" 1 ) 
and again shows only slight variation between the models, 
with a FWHM that is now dominated by the thermal width 
because of the lower mass of H. 

The variation in the properties of the synthesized line pro- 
files as the models evolve is presented in Figure ^3 which 
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Fig. 10. — Evolution of emission line kinematic properties for a face-on 
viewing angle for the three models. The different models are shown by: A 
(solid line), B (dashed line), and C (dotted line), (a) Mean line velocity, v. 
(b) Root-mean-squared line width, <x. In panel (a), the top set of lines are 
for the neutral collisional line, the middle gray set are for the recombination 
line, and the bottom set are for the ionized collisional line. In the panel (b) 
the order is reversed. 



shows the mean velocity, 



Jo 

v = ^ 



(z)t](z) dz 



(7) 



rj(z) dz 



and the RMS velocity width, 



(v(z)-v 

Jo 



fr](z)dz 



(8) 



rj(z) dz 



calculated for the three models as a function of time. Note 
that the velocity width, cr, only includes the contribution from 
the acceleration of the gas and does not include the thermal 
Doppler broadening. Supposing the emission to come pre- 
dominantly from gas at 10 4 K, the observed FWHM should 
be related to the RMS width by FWHM [(20 km s" 1 /A) 2 + 
(2.3550-) 2 ] 1 / 2 . 

Once the models have settled down to their quasi-steady 
stage of evolution, it is remarkable how little variation is seen 
in the emission line properties, either between the models or 
as a function of time. 

4. OBSERVATIONAL PARAMETERS OF THE ORION NEBULA 

In this section, we derive observational parameters for the 
Orion nebula in order to compare with our simulations. We 
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choose to use the spatial variation of the emission measure 
and electron density, together with the kinematic properties 
of various optical emission lines. 

4.1. Emission measure and electron density 

Emission measure and electron density maps of the inner 
Orion nebula are shown in Figure[^ We base these maps on 
velocity-resolved dat a cubes of the emissivity in various op- 
tical emission lines ( Garcia Diaz & Hennevll2003t iDoi et all 
2004). This allows us to distinguish between the emission 
from the photoe vaporation flow an d that produced by Herbig- 
Haro shocks JO' Dell et al.lll997h and back-scattering from 
dust in the molecular cloud (Hennev 1998). 

The emiss i on m easure map is calculated from the Ha data 
of D oi etal.1 (|2004), using a heliocentric velocity interval of 
+ 1 to +33kms _1 , and correct ing for foreground dust extinc- 
tion using the C(Hf3) map of O' Dell & Yusef-Zadehl (2000) 
and the extinction curve of Cardelli et al. ( 1989). The emis- 
sion measure is assumed to be strictly proportional to the Ha 
surface brightness, which would be the case for isothermal, 
Case B emission, and the constant of proportionality is found 
by co mparison with the measured fr ee-free optical depth at 
20 cm (lO'Dell & Yus ef-Zadeh 2000), assuming an electron 
temperature of 8900 K. In regions of the map that do not con- 
tain significant high-velocity emission our derived emission 
measure is identical to that derived directly from the radio 
free-free data. 

As an additional check on our results, we flux-calibrated 
our Ha map (using the full velocity range ) by comparison 
with published H ST WFPC2 imaging do' Dell & DoilH999l 
iBallv et al1l2000l) and co nverted to emission measure using 
the atomic parameters in Osterbrock ( 1989), which yielded 
value s consistent at the 10% level with those derived by 
iBaldwin etal.1 dl991 | ) from the H 1 1-3 line. As recognised 
previously dWen & Q'DelilfT 995). the emission measures de- 
rived from the full velocity range of optical recombination 
lines are overestimated by up to 50% because of the contri- 
bution of light that has been back-scattered by dust in the 
molecular cloud. Our velocity-resolved technique automati- 
cally compensates for this without the need to introduce an 
arbitrary correction factor. 

In order to calculate the electron density in the pho- 
toevaporation flow, we use velocity-resolved maps of 
the density-sensitive doublet lines [S ii] 6716 and 6731 A 
( iGarcfa Diaz & Hennevll2003L 120051) to select only that gas 
in the range +4 to +24 km s _1 , close to the peak velocities of 
fully ionized species. The summed surface brightness in this 
velocity range for the two lines is shown in Figure [TTt> after 
correction for foreground extinction (the apparent emission 
peak at the lower right of the image is an artifact of the reduc- 
tion procedure). The density-sensitive ratio R = 6731/6716 
for the same velocity range is shown in Figure Hit, in which 
higher values of R correspond to higher densities. In order to 
calibrate our ratios, we first calculated maps summed over the 
full velocity range of the [S n] emissio n and compared t hese 
with t he spectrop hotometric data of Baldwin et al.] ( 119911) and 
Pogg e et al1dl992l) . 

We find that a very good fit to the dep endence of electron 
density on R for gas at 8900 K dCai & Pradhanl 19931) is 



where n = 2489cm- 3 , Ri = 0.697, and R 2 = 2.338. Us- 
ing this fit, the smallest ratio shown in Figure II lb corre- 



sponds to n e ^ 1000 cm 3 and the mid-range corresponds 
to 4-6000 cm -3 . The largest ratios shown in the image are 
in the high-density limit (n e > 10 4 cm~ 3 ), at which point the 
derived density becomes very uncertain due to the extreme 
sensitivity to experimental errors and variations in the tem- 
perature. Previous work, which did not resolve the kinematic 
line profiles, fo und that electron densities derived from fully 
ionized species dJonesll 1992L IMegeath et al.lH990l) are gener- 
ally higher than those determined from [S n]. This has been 
attributed to the fact that part of the [S n] emission comes from 
partially ionized zones in which the electron fraction is signif- 
icantly less than unity. As a result, it was found necessary to 
apply a correction factor to the [S n] densities, although there 
has been significant d isagreement on its value (Baldwi n et alJ 
fl99Tt IWen & O'DelH^ . By using only the blue-shifted 
portion of the line profiles to calculate our densities, we are 
selecting the [S n] emission from the fully ionized gas. As 
a result, we have no need to apply a correction factor to our 
densities. Indeed we find that the 6731/6716 ratio for the 
blue-shifted velocity range is typically 10% larger than the 
ratio calculated for the entire line, although there are some 
local variations and the difference is less at greater projected 
distances from the ionizing star. 

4.2. Gas kinematics 

In order to compare with the gas-kinematic behavior of the 
numerical models, we here present kinematic maps of the op- 
tical emission lines [O i] 6300 A and [S m] 63 12 A (FigurefHl. 
These are derived from echelle spectroscopy through a NS- 
oriented slit at a series of 60 different pointings using the Mez- 
cal spectrograph of the Observatorio Astronomico Nacional, 
San Pedro Maitir, Mexic o. The observation s are described 
in more detail in Garcia Diaz & HennevH2005l) . We use these 
data in preference to the [N n], Ha, and [O in] data of D oi et alJ 
(2004) because the latter contain no information on the east- 
west variation of the line velocities since they assume that the 
peak velocity of the total profile summed over each NS slit is 
the same. With our data, on the other hand, we have taken 
a pair of EW-oriented slit spectra, which allow us to "tie to- 
gether" the wavelength calibrations of the individual NS slits. 
In addition, the existence in our spectra of the night-sky com- 
ponent to the [O i] line (at rest in the geocentric frame) allows 
us to make a very accurate wavelength calibration for that line. 
The images in Figure[2]show the spatial variation in line sur- 
face brightness (which has not been corrected for foreground 
extinction), the intensity-weighted mean line velocity, and 
the non-thermal component of the full-width half-maximum 
(FWHM) of the line, which was calculated by the subtrac- 
tion in quadrature from the observed width of the instrumen- 
tal width (9kms~') and the thermal width (5.35 km s~' for 
[Oi], 3.78 km s -1 for [S m]). All quantities were calculated 
for a spectral window that extends between heliocentric ve- 
locities of -12 to +40kms _1 to minimize the contamination 
from high-velocity emission from shocks. 8 

The [Oi] line i s an example of a neutral collisional line, 
as discussed in § 13.41 As such, it originates in a thin zone 
close to the ionization front and shows a lot of fine-scale struc- 
ture in its surface brightness. The mean velocity of [Oi] is 
only blueshifted by a few km s _1 with respect to the molecu- 

8 Although this is successful at eliminating the highly blue-shifted emis- 
sion associated with Herbig-Haro jets, there is also intermediate-velocity 
emission associated with the bowshocks of these objects that does enter our 
velocity window and skews the velocity maps in localized patches. 
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lar gas, which emi ts at 25-28 km s 1 ( Omodaka et alJI 19941 
Wils on et al.ll200"l . There is a general trend for the [Oi] 
emission to become more blue-shifted towards the west of 
the Trapezium, which is also seen in CO emission from the 
molecular gas (Wilson et al. 2001). Some of the bright lin- 
ear filaments seen in the surface brightness map (such as the 
Bright Bar in the SW) correspond to filaments of increased 
redshift in the mean velocity map. Other linear structures 
in the mean vel ocity have no visible co unterpart in the sur- 
face brightness (Garcia Diaz & Hennev 2005). The most red- 
shifted [O i] comes from diffuse emission just above the SE 
end of the Bright Bar, while the most blueshifted is due to 
the residual effect of high-velocity shocked emission from 
HH201, 202, 203, and 204. The [Oi] line width shows little 
large scale structure, with the notable exception of a marked 
decrease in the line width along the Bright Bar filament. Lo- 
calized spots of increased line width correspond to the HH 
objects discussed above. 

For the ionization parameter and EUV spectrum relevant 
to Orion, the [S m] emission traces the fully ionized gas 
(x > 0.9), except for the zone closest to the ionizing star, 
which is predominantly S 3+ . It comes from a thicker layer 
and so shows smoother variations in surface brightness than 
does [Oi], being similar in appearance to the emission mea- 
sure map of Figure^J The mean line velocity is blue-shifted 
several kms 1 with respect to [Oi] and shows a more pro- 
nounced gradient from east to west. On the other hand, the 
large-scale features of the mean velocity map are rather sim- 
ilar to [Oi], with redder emission being associated with the 
emission from the Bright Bar and other linear features. One 
feature seen in [S m] but not in [O i] is an EW-oriented spur 
of blue-shifted emission to the north of the Bright Bar (the 



O m counterpart of this feature is discussed in lDoi et alJ2004l) . 
The average [S m] linewidth is not much different from that of 
[O i] but [S m] shows much more variation across the nebula. 
The highest widths (15-20 km s -1 ) are seen near the Trapez- 
ium, at the east end of the blue spur, and in the red-shifted 
emission at the faint SW end of the Bright Bar. The lowest 
widths (^ 12kms~') are seen along the length of the Bright 
Bar and in a ring of radius 30-60" around the Trapezium. 

5. APPLICATION OF THE MODELS TO ORION 

In order to compare our numerical simulations with the ob- 
servational data of the previous section, we first extract the 
observational data in an EW strip, centered on the ionizing 
star 9 l Ori C. The strip is 60" wide in the NS direction 
and the average emission measure (calculated from the Ha 
surface brightness) and electron density (calculated from the 
[S n] line ratios) are shown in Figure^]as a function of pro- 
jected distance from the ionizing star (assuming a distance of 
43Qpc. lO'Dell (l2001bl) V The data are represented by a gray 
band that indicates +3cr variation of individual pixels from 
the mean. The derived electron density is very uncertain close 
to 6 l Ori C on the West side since the observed line ratio is 
close to the high-density limit. We have therefore omitted this 
data from the graph. 

We use the three models A, B, and C at a common age of 
78,000yr as representative of our numerical results and cal- 
culate the simulated observed profiles of emission measure 
and electron density for a viewing angle of 15° with respect 
to the line of sight, following a procedure as close as possible 
to that employed in calculating the observed maps above. For 
the simulated emission measure, we first calculate simulated 
Ha maps, restricted to the line-of-sight velocity range -4 to 
-23 km s _1 , then convert this to emission measure by assum- 
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Fig. 12. — Emission line maps for [O I] 6300 A (top row) and [S III] 63 1 2 A (bottom row) constructed from longslit echelle spectroscopy. Left panels show the 
surface brightness. Central panels show the mean heliocentric line velocity. Right panels show the FWHM line width corrected for thermal and instrumental 
broadening. 
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Fig. 13. — Observed East-West profiles of emission measure (top panels) 
and electron density (bottom panels) in the Orion nebula (gray line), com- 
pared with predictions of the model simulations for an inclination angle of 
15°: Model A (solid line), B (dashed line), and C (dotted line). 



ing a fixed Ho- emissivity corresponding to 8900 K. For the 
electron density, we calculate simulated maps for the 6716 A 
a nd 6 731 A [Sn] doublet components using the methods of 
§ !3.4l for the same line-of-sight velocity range as for Ha. We 
then use the ratio of these two maps to derive « e in exactly the 
same way as for the observations. 9 

For each model, we rescale the lengths and densities to the 
Orion values as follows: first, we determine the star-ionization 
front distance, Zh,o> by comparing the FWHM of the width 
of the EM in model units with the observed value of 6.12 x 
10 cm. Then, we determine the electron density scaling by 
fitting by eye to the broad peak of the observed EM profile. 

The results are shown in the upper panel of Figure ^] in 
which it can be seen that all three models can provide a rea- 
sonable fit to the emission measure profile for a viewing an- 
gle of 15°. Models A and B fit the large-scale skewness of 
the profile somewhat better than Model C, which overpre- 
dicts the EM on the West side and underpredicts it on the East 
side. None of the models can reproduce the fine-scale struc- 
ture seen in the observed profile. In general, the agreement is 
better on the West side of the profile, which is smoother than 
the East side. Model profiles for a viewing angle of 30° (not 
shown) badly fail to fit the observed profile since their peaks 
are shifted too far to the West, while those for a viewing angle 
of 0° (not shown) do not fit as well as the 15° profiles since 

9 Since the ionization fraction of S + does not exactly follow that of either 
neutral or ionized Hydrogen, we have used a simple fit to the dependence of 
of the S + on x that we derive from static photoionization equilibrium models 
calculated using the Cloudy plasma code iFerland 2000): 

«(S + ) _ q(l - xf 
n(S) ~ 1 + (a - 1)(1 - xf ' 

with a = 6.16, b = 0.730, c = 0.629. 



TABLE 1 

Model fit parameters for 15° viewing angle 



(WMh/IO^s- 1 Zh , /10 17 cm Zp ,o/10 17 cm n p /10 4 c m - 3 



Model 


(1) 


(2) 


(3) 


(4) 


A 


0.299 


4.42 


3.80 


0.707 


B 


0.363 


3.88 


3.55 


1.108 


C 


0.919 


4.93 


4.65 


1.411 



(1) Effective stellar ionizing luminosity. (2) Distance from star to heating 
front. (3) Distance from star to density peak. (4) Peak ionized density. 
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Fig. 14. — Observed East- West profiles of [O I] (mid gray band) and [S III] 
(light gray band) mean velocity (top panels) and FWHM (bottom panels) 
in the Orion nebula, compared with predictions of the model simulations: 
Model A (solid line), B (dashed line), and C (dotted line). 



they are too symmetrical. 

The parameters for the model fits in the 15° case are listed 
in Table [2 These are the required effective ionizing photon 
luminosity of the central star, the distances along the symme- 
try axis from the star to the heating front and to the peak in 
the ionized gas density (see §|2ji, and the value of the ionized 
density peak. Since dust grains in the ionized nebula will ab- 
sorb a f raction, fd ~ 0.5, of the ionizing photons emitted by 
the star jArthur et alJl2004l) . and this effect is not included in 
the simulations, the derived quantity is (1 - fdjQn- 

The results for the simulated observed electron density are 
shown in the lower panel of Figure fT3*l where it can be seen 
that Model C does the best job of fitting the observations. The 
other two models significantly underpredict the density in the 
central and Western portions of the nebula. None of the mod- 
els reproduce the very high (but poorly determined) de nsities 
that seem to be indicated to exist close to 8 1 Ori C. In § 16.1.31 
below, we discuss how the the geometry of Model C can be 
reconciled with the fact that the nebula appeal's concave on 
the largest scales. 
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We have also calculated simulated maps of mean velocity 
and line width of the [O i] and [S m] lines in order to compare 
with the observed kinematics presented in § 14.21 The results 
are shown in Figure^Jfor the same EW strip as was used for 
the emission measure and density comparisons. The upper 
panel shows the observed mean heliocentric velocities of the 
two lines as gray bands indicating +3cr variation of individual 
pixels at each RA. The lines show the model predictions for 
an inclination angle of 15° and assuming that the molecular 
gas has a constant heliocentric velocity of +28kms~ 1 . The 
lower panel shows the same but for the non-thermal FWHM 
of the lines. 

The models reproduce the magnitude and trend of the [O i] 
velocity as well as the qualitative fact that the [S m] is more 
blue-shifted than [Oi]. Although they agree with the [S m] 
mean velocity on the western side of the nebula, they some- 
what overpredict the blueshift of the [S m] line in the center 
of the nebula. Overall, Model C is a closer fit than the other 
models. 

With regard to the line widths, the agreement is not very 
good. The observed [S m] widths are 2-5 km s -1 larger than 
the best model predictions (from Model C) and the [Oi] 
widths are 7-9kms~ 1 larger than any of the model predic- 
tions. 

6. DISCUSSION 
6.1. The Orion Nebula 

In § |3 the model predictions for the profiles of emission 
measure, electron density, and line kinematic properties are 
compared with observations of the Orion nebula. The EM 
profile is used to calibrate the length and density scales of the 
models and so its comparison with the model profiles does 
not provide a strong test. It is worth noting, however, that 
Models A and B can reproduce the large-scale skewness of 
the EM profile for a viewing angle of 15°, whereas Model C 
cannot. For Model C to produce the observed skewness re- 
quires a larger viewing angle, which shifts the predicted EM 
peak much farther to the West of 9 l Ori C than is observed. 
A further discrepancy with all the models is that the observed 
profile shows much fine-scale structure that is entirely absent 
in the simulations. 

Since the observed n e profile was not used in determining 
the model parameters, it provides a much more stringent test 
for the models. Only Model C predicts a sufficently high den- 
sity in the central and Western regions of the nebula. Again, 
the observations show structure on scales < 10 17 cm, which 
are not present in the models. 

6.1.1. Required ionizing photon luminosity 

The derived parameters from the model fits (Tabled show 
striking differences between the three models. In particular, 
the required effective ionizing photon luminosity is 2h(1 - 
/ d ) = 9.2 x 10 48 s _1 for Model C (where / d is the fraction of 
ionizing photons absorbed by dust), which is roughly three 
times higher than the equivalent values for Models A and B. 
This is partly because the ionized density is more concen- 
trated at the ionization front in Model C and thus requires 
a higher ionizing luminosity to sustain the same emission 
measure as compared to a situation in which the ionized gas 
is more evenly distributed between the star and the ioniza- 
tion front, as is the case for Models A and B. Another con- 
tributing factor is that the convexity of the ionization front 
of Model C (see Figure [3J means that the nebula is density- 
bounded for a large fraction of the solid angle as seen from 



the central star, whereas Models A and B are more concave 
and thus more ionization bounded. It has long been known 
that neutral gas in the so-called "veil" overlays the Orion neb- 
ula ( Ivan der Werf & Gosslll989tlAbel et"dE0 04). so that the 
nebula must be ionization bounded on a scale of a few parsecs, 
thus contributing a low surface-brightness halo to the ionized 
emission. However, this would make a negligible contribu- 
tion to the brightness of the inner region studied here. Indeed, 
Felli et al. ( 1993) find the total effective ionizing luminosity 
driving the nebula to be 9 x 10 48 s _I from single-dish radio 
observations, which are sensitive to emission on a scale of par- 
secs that is not detected by interferometers. This is fully con- 
sistent with Model C but not with the other two models. An- 
other independent estimate of the stellar ionizing luminosity 
comes from stud ying the surface brightnes s of the cusps of the 
Orion proplyds llHennev & Arthur! 1 19981 iHennev & O'Defi 
I1999I) . which gives a value of (1- f'^On = (1.0±0.2)xl0 49 s _1 
for proply d 167-317 (LV2) th at has the best-determined pa- 
rameters llHennev et alJ |2Q02). The fraction of ionizing pho- 
tons absorbed by dust in the proplyd flow, ft, should be sim- 
ilar that for the nebula as a whole, /<j, since in both cases the 
dust optical depth through the ionized gas is of order unity, so 
this result is also consistent only with our Model C. 

6.1.2. Emission line kinematics and broadening 

The comparison between the predicted and observed line 
kinematics is only partially successful. On the one hand, 
the trend of increasing blueshift with ionization is well re- 
produced by the models. This is a vast impr ovement on 
the pla ne-parallel weak-D models presented in Hennev et alJ 
(2005), which were unable to produce a difference of more 
than lkms -1 between the neutral and fully ionized lines. 
Our photoevaporation models predict a difference of about 
5kms~' between [Oi] and [Sm], which is very similar to 
what is observed. The observations show a shift in the mean 
velocity from E to W of ^ 3kms _1 for [Oi] and =s 6kms~ 1 
for [Sm], which are both approximately reproduced by the 
models. In the case of [Oi], some part of this variation can 
be accounted for by shifts in t he velocity of the m olecular 
gas behind the ionization front ( Wilso net alJl200il) and fur- 
thermore the westernmost portion of the [S m] is affected by 
intermediate-velocity gas associated with the HH 202 jet and 
bowshock. Taking these into account, Model C can be seen 
to fit the observations quite well, apart from in the region 
around the Trapezium where the observed [S m] velocity is 
about 4 km s _I redder than the model prediction. This may be 
evidence for the interaction of the stellar wind from 9 1 Ori C 
with the photoevaporation flow, as is discussed further below. 

On the other hand, the discrepancy seen between the pre- 
dicted and observed line widths is very striking, especially for 
[Oi], for which the models all predict a FWHM of 7kms _1 , 
which is only half the observed value. For [S m], the discrep- 
ancy, although less, is still significant: the best model predicts 
a FWHM of 12kms~' whereas the observed values range 
from 1 3 to 17 km s~ 1 . Since the widths of independent broad- 
ening agents should add approximately in quadrature, we thus 
require an extra broadening of 12 km s _I for [O i] and 9 km s 
for [Sm]. The line widths in the models are mainly deter- 
mined by the physical velocity gradients, dv/dz, in the pho- 
toevaporation flow and are of magnitude ~ \dv/dz\Sz, where 
Sz is the thickness of the emitting layer. In principle, an ad- 
ditional broadening agent could be variation along the line of 
sight of the direction of the flow velocity vector (with magni- 
tude ~ |y| 89) but this effect is very small in our models since 
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the radius of curvature of the ionization front is of order the 
size of the nebula, so 50 is small. This broadening would be 
much more important if the ionization front had many small- 
scale irregularities, which does seem to be indicated by the 
large amount of structure seen in the surface b rightness maps 
( Hester et al] ll991b lO'Dell et alJl!9911 120031) . An example 
of where this effect has been sho wn to be important is in 
the flows from the Orion proplyds jHennev & Q'Dell fl999). 
These are poorly resolved spatially with ground-based spec- 
troscopy, so the full range of angles is present and one predicts 
line widths of order 2 |d|, with photoevaporation models suc- 
cessfully reproducing the observed line profiles for fully ion- 
ized species. However, although this small-scale divergence 
broadening might plausibly explain the [S m] linewidth dis- 
crepancy, it is hard to see how it can work for [O i] since the 
magnitude of the velocity in the [O i] -emitting zone is only 2- 
3 km s _1 , which, even when doubled, is much smaller than the 
required broadening of 12 km s -1 . 10 Similar arguments apply 
to the contribution to the observed widths o f back-sca ttered 
emission from dust behind the emitting layer tHennevI 19981) . 
which again may be important for blueshifted lines such as 
[S m] but will be unimportant for [O i]. 

Previous authors have reached similar conclusions about 
the need for an extra broadeni ng agent in the nebula from a 
purely empirical viewpoint (see lO'De l]|2001b [ 8 3.4.2- 4 and 
references therein). More recently, O'Dell et al. (2003) have 
found evidence that the extra broadening may be greater in 
recombination lines than in collisional lines and thus may be 
related to the equally puzzling problem of temperature varia- 
tions in the nebula. 

Another proposed solution is that the broadening is caused 
by Alfven w aves, as is believed to be th e case in giant molec- 
ular clouds (M vers & C.oodmanl [T988). In order to give an 

Alfven velocity (u a = (B 2 l4np} ^ ) of 9kms _1 , which is nec- 
essary to explain the observed [O i] broadening, one requires 
a magnetic field of order B = 0.5 mG in the region where the 
line forms (which has an ionization fraction of order 20% and 
a total gas density roughly 3 times the peak ionized density). 
This is marginally inconsistent with the upp er limit on B de- 
rived from the Faraday rotation measure ( Ra o et alJ 19981) and, 
furthermore, would predict very high Alfven speeds either in 
the dense molecular gas or in the lower density, highly ionized 
gas. To see this, consider two possibilities for how the mag- 
netic field responds to compression/rarefaction: first, suppose 
that B is constant, such as might be the case for a highly or- 
dered field with orientation parallel to the flow direction (per- 
pendicular to the ionization front). In this case, va ~ p~ 1 ^ 2 and 
so a reasonable value is obtained in the dense molecular gas 
(va - 3 km s _I at 10 5 cm -3 ) but a very high value is obtained 
in the more highly ionized, lower density regions of the flow 
(va — 30km s at 10 3 cm -3 ), which is ruled out by observa- 
tions of lines such as [O m]. Second, suppose that B is tangled 
or chaotic and can be considered as an is otropic pressure with 
effective adiabatic index y m = 4/3 (see Hennev et al J 120051 
Appendix C). In this case, v A ~ p (rm ~ I)/2 ~ p 1/6 , so that the 
Alfven speed would have to exceed lOkms^ 1 in the molecu- 
lar gas, which greatly exceeds the linewidths of the molecu- 
lar lines. Although we cannot definitely rule out the role of 
Alfven waves, these considerations suggest that they do not 
provide a natural explanation for the line broadening. 

10 The [Ol] linewidths from the proplyds are also inexplicable in terms 
of kinematic broadening since t hey are very similar to the w idths in the sur- 
rounding nebula I H ennev & O'D ell 1999; O'Dell et alJ200!l) . 



6.1.3. Consistency with the large-scale nebular geometry 

Comparison of our models with the observed electron den- 
sities and line velocities in the nebula strongly favor Model C 
over the other two models, so it is worth considering whether 
the shape of the ionization front in this model is consistent 
with independent determinations. In Model C, the ionization 
front has a positive curvature, that is, it is convex as seen from 
the ionizing star. However, three-dimensional reconstruction 
of the shape of t he ionization front paints a somewha t more 
complex picture (W en & OT5elilll995l iHennevI 12005). as is 
shown in Figure^] In the brightest part of the nebula, just 
to the W of the ionizing star, the ionization front is indeed 
locally convex in the EW direction but is roughly flat in the 
NS direction. 1 1 Also, the axis between the star and the closest 
point on the front is inclined at about 15° to the line of sight, 
consistent with what we found for our models. On a larger 
scale, however, the front becomes concave in all directions 
apart from the south-west. The implications for our models 
are twofold: first, on a medium scale (< V ^ 4 X 10 17 cm) 
the flow should not be strictly cylindrically symmetric but in- 
stead should resemble Model C in its EW cross-section but 
more resemble Model B in its NS cross-section. Since we 
compared our models against an EW strip of the nebula, it is 
reasonable that Model C produced better fits. Secondly, on a 
larger scale the flow should come to more resemble Model A, 
but with a symmetry axis that is inclined towards the south- 
west, rather than towards the east and possible at a greater 
angle to the line-of-sight. 12 How the flow will adjust between 
the differing curvature and symmetry axes is unclear, but it 
may involve internal shocks. Further complications will arise 
due to the interaction of the global flow with local flows from 
other convex features in the ionization front. The largest and 
most prominent of these is the Bright Bar but there are many 
other examples (O'Dell & Yusef-Zadeh 2000), some of which 
are apparent in Figure [21 

6.2. Shortcomings of our models 

Although our photoevaporation flow models have had some 
measure of success in accounting for the structure and kine- 
matics of the Orion nebula on medium scales, there are many 
physical processes that we have neglected and which would 
need to be included in a complete model. 

We have used a very simplified prescription for the radia- 
tive transfer in our models, in which we consider only one 
frequency for the ionizing radiation. Although we do include 
the hardening of the radiation field, thus accounting for the 
important reduction in the photoionization cross section deep 
in the ionization front, it would be more satisfactory to use 
a multi-frequency approach. The consideration of higher en- 
ergy EUV radiation would allow us to calculate the helium 
ionization structure, while the inclusion of non-ionizing FUV 
radiation would permit a more realistic treatment of the heat- 
ing and cooling in the neutral gas. Furthermore, the effects of 
dust absorption on the transfer of ionizing radiation are im- 
portant and should be included. 

One failure of our model has been an inability to acount 
for the reduced blueshift of high-ionization lines that is seen 
within a few times 10 17 cm of the Trapezium. It is possible 

" The front seems to wrap around the dense molecular filament that is 
seen in molecular line and dust continuum observations j Johnston e & Balhl 
1999). 

12 This is also consistent with the parsec-scale optical appearance of the 
nebula, in which the flow appears to be towards the south-east. 
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Fig. 15. — Three-dimensional shape of the principal ionization front in Orion, derived from radio free-free surface brightness i Hennev 2005). (a) View over the 
molecular ridge from the WNW direction, (b) View over the Bright bar from the SSE direction. The x, y, and z axes correspond to displacements along the EW, 
NS, and line-of-sight directions, respectively, with numerical values marked in arcseconds (1" = 6.45 X lO 15 cm). The ionizing star is at the origin and is marked 
by a cross, while the observer is at (x, y, z) = (0, 0, +oo). Shading of the surface indicates the z height from low (red) to high (blue). This figure is available in 
color in the online edition. 



that the stellar wind from 6 l Ori C or radiation pressure acting 
on dust grains may be important in this inner region. Both are 
effects that need to be incorporated in a future model. 

There is muc h small scale struct ure present in the Orion 
ionization front jHesteretalJ ll991). right down to scales ~ 
10 15 cm, 13 or smaller still if one includes the proplyds. Our 
models can be taken as reflecting the mean flow properties 
after all this small scale structure has been smeared out, but 
it is not obvious that such an approach is strictly valid. The 
very highest electron densities found close to the Trapezium 
are not reproduced by our models and seem to require that the 
emission there is dominated by much smaller-scale flows than 
we have considered. The effects of the mutual interac tion of 
such flows is still in need of further research (but see lHennevl 
2003), which also needs to address the influence of bar-like 
features and the change in the orientation of the symmetry 
axis when one passes from 10 18 cm to 10 19 cm scales. 

Lastly, but perhaps most importantly, our models com- 
pletely neglect the role of the magnetic field. In the cold 
molecular gas, the pressure associated with magnetic field is 
expected to dominate the gas pressure, and so should have an 
important dynamic effect, but it is unclear if this is still the 
case in the ionized gas. For most plausible assumptions about 
the magnetic field geometry, the magnetic pressure is found to 
be a small fraction o f the thermal pressure in the ionized gas 
dHennev et alJl2005l) . If the unexplained broadening of the 
[O i] line found in § |4.2l is ascribed to Alfven turbulence, then 
this fraction may be as high as 25% at the ionization front. 
However, even if the magnetic field is dynamically unimpor- 
tant in the ionized gas, it can still have an important effect 
on the jump conditions at the ionization fron t ( Redma n et al.l 
H998tlWilliams et alJ2000HWiUiams & Dvsonl2001l) . " 

7. CONCLUSIONS 



We have carried out two-dimensional, cylindrically sym- 
metric, hydrodynamical simulations of the photoevaporation 
of a cloud with large-scale density gradients, which gives rise 
to a transonic, ionized, photoevaporation flow with the fol- 
lowing general properties: 

1. After an initial transient phase, with dur ation ~ 10 4 yr 
for typical compact H n region sizes (§ 13. 1> . the flow 
enters a long-lived quasi-stationary phase in which the 
ionized flow is approximately steady in t he f rame of 
reference of the ionization/heating front fS 13.21 , 



During the quasi-stationary phase (§§ 13. 2113731 . the flow 
structure is determined entirely by two parameters: the 
distance of the ionizing star from the front, zo, and 
the curvature of the front, k. Depending on the cur- 
vature, the flows are more "champagne-like" (a: < 0) or 
"globule-like" (k > 0). 

The curvature of the front (§ 13. 3> depends on the lat- 
eral density distribution in the neutral gas and on the 
evolutionary stage of the flow. Positively curved fronts 
are only obtained for lateral density distributions that 
are asymptotically steeper than r~ 2 . We found that ini- 
tially flat or concave fronts tend to evolve with time 
towards more negative values of kzo (increasing con- 
cavity), whereas initially convex fronts evolve towards 
more positive values of kzq (increasing convexity). 

The flat and convex ionization fronts in the simulations 
are found to be D-critical, whereas concave fronts are 
found to be weak-D (AppendixlBl. There is still a sonic 
point in the flows from concave fronts but it occurs in 
fully ionized gas, rather than in the front itself. 



13 It is unclear whether this structure reflects pre-existing structure in the 
molecular cloud or whether it is due to instabilities associated with the ion- 
ization front I Garcia-Segura & Franco 1996; Williams 1999, 2 003) . 



5. The relationship between gas ionization and kinematics 
sh ows little variation as a function of the front curvature 
(§ I3.4> . In all cases, a gradient of lOkms -1 between 
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the neutral and ionized gas is obtained. There is, how- 
ever, a tendency for models with more positive curva- 
ture to show larger mean velocities and line widths for 
ionized species. 

We have compared our models with both new and existing 
observations of the Orion nebula (§§EI|5), with the following 
results: 

1 . Only a convex model can simultaneously reproduce the 
observed distributions of emission measure and elec- 
tron density in the nebula. A good fit is obtained for 
axo — 0.5, implying that the radius of curvature of the 
front is roughly twice its distance from the ionizing 
star. Flat and concave models predict an electron den- 
sity in the west of the nebula that is much lower than 
is observed. They do, however, fit the observed emis- 
sion measure distribution slightly better than a convex 
model. 

2. Only a convex model is consistent with independent es- 
timates of the ionizing luminosity, Qu, from Q l Ori C. 
Flat and concave models trap a much higher fraction of 
the stellar ionizing photons in the central region of the 
nebula and predict a Qu that is three times too small. 

3. All the models can broadly reproduce the observed 
mean velocities of optical [O i] and [S m] emission lines 
but a convex model fits the data better. 



4. None of the models can reproduce the observed widths 
of the optical emission lines, although convex models 
come close to doing so for the high-ionization lines. 
The largest discrepancy is found for the neutral [O i] 
line, which requires an extra broadening agent with a 
FWHMof 12kms'. 



5. The best-fitting model over all has a flow axis inclined 
15° west from the line of sight and a distance between 
8 l Ori C and the ionization front of 4.7 x 10 17 cm. 
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Fig. A16. — Schematic diagram of the photoevaporation flow from a plane ionization front illuminated by a point source, in which the principal geometrical 
quantities are indicated. 
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APPENDIX 

A. ANALYTIC MODEL FOR THE FLOW FROM A FLAT IONIZATION FRONT 

In this section, we develop an approximate analytic treatment of the ionized portion of the photoevapor ation flow, for the 
special case in which the ionization front is exactly flat. An earlier version of this calculation was presented in lHenne 



A.l. Assumptions and Simplifications 

Consider an infinite plane ionization front (z - 0) illuminated by a point source of ioniz ing radiation, which lies at a height 
z — z* above the ionization front and defines the cylindrical axis r — (see Figure IaToY Suppose that the ionization front 
is everywhere D-critical and that the volume z > is filled with a time-steady isothermal ionized photoevaporation flow with 
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straight streamlines. Each streamline can be labeled with the cylindrical radius ro of its footpoint on the ionization front and 
makes an angle /(ro) with the z-direction. Assume that the acceleration of the flow along a streamline is identical to that in the 
photoevaporation flow from a spherical globule with a radius equal to the local "divergence radius" of the flow: R d = ro esc / so 
that the velocity along each streamline follows an identical function U(y) = u/ci, where y — 1 + s/R d , s is the distance along the 
streamline, and c; is the isothermal sound speed in the ionized gas (since the front is D-critical, we have t/(0) = 1). The density 
at any point (r, z) can therefore be separated as 

N(rp) 

w(r ' z) = 2777 n (A1) 

where r = ro + (y - l)Rd sin i, Z = (y - l)R d cos U an d N(ro) is the density at the ionization front. 

The effective recombination thickness, h, is defined as in the caption to Figure except integrated along a streamline rather 
than along the z-axis: N 2 (ro)h = j n 2 (ro, s) ds, yielding 



h = R d j y^U^dy = a)R d , (A2) 

where a> is a constant that depends on ly on the acceleration law U(y). Although a> 0.1 is typical for spherical photoevaporation 
flows from globules (Bertoldil ll989l) . the results of our numerical models above (see Fig.0 suggest that the acceleration is 
somewhat slower in the present circumstances, leading to a larger value of h, so we leave u> as a free parameter for the time being. 

A. 2. Photoionization Balance 

With this definition of h, and assuming that the photoevaporation flow is "recombination dominated" llHennevl UoOl). that 
recombinations can be treated in the on-the-spot approximation, and ignoring any absorption of ionizing photons by grains, one 
can approximate the recombination/ionization balance along a line from the ionizing source to the ionization front by considering 
the perpendicular flux of ionizing photons incident on the equivalent homogeneous layer: 

Oh cos 6 1 

- a B N\r )h(ro), (A3) 



4n(zi + rQ 

where tan# = ro/z*, Qh is the ionizing photon luminosity (s _l ) of the source, and is the Case B recombination coefficient. 
This approximation will be valid when h is small enough to satisfy h < z* and h sec 6 < H, where H = \dltiN/dro\~ l is the radial 
scale length of the ionized density profile at the ionization front. This profile can be found from equation jA3t to be 



N(r ) = N 



2\ 



>4 



- 3/4 I i \"l/2 

' n 1 



(A4) 



where x = ro/z, and A^o, ho are the values on the axis at ro = 0. Note that for any plausible h(ro), the density at the ionization 
front will decrease with increasing radius. 

A.3. Lateral Acceleration 

In order to complete the solution, it is now only necessary to specify the radial dependence of the streamline angle, /(ro), 
from which the radius of divergence, Rd(ro) and scale height, h(ro), follow automatically. To do this, we must recognise that 
the streamlines cannot truly be straight: in a D-critical front the gas will be accelerated up to the ionized sound speed, c;, in the 
z-direction (perpendicular to the front) but at z = it will have no velocity in the r-direction (parallel to the front). However, the 
pressure gradient associated with the radially decreasing density will laterally accelerate the gas in the positive r-direction as it 
rises through the recombination layer (which is much thicker than the front itself), thus bending the streamlines outwards. To be 
definite, we will consider the characteristic angle of each streamline, /(ro), to be that obtained by the gas at z = h( r o)- 

The lateral acceleration, a, of the gas is given by 

1 dP c 2 

a = if£ = _L ; ( A5) 
p dr H 

where H is the lateral density scale length defined above. Immediately after passing through the heating front, the gas will have 
a vertical velocity of 0.3ci, rising to c; once it is fully ionized and then to 2q as it passes through the recombination 
layer. We take the typical vertical speed to be v z - Bci, where B is a parameter whose value must be determined but which 
should be somewhat less than unity. Hence, the time taken to rise to a height h is St — h/Bc[, so that the lateral speed reached is 
u r = adt = cJijBH and the streamline angle is given by 

v r h 

tan / = — = — , (A6) 
v z BH 

so that R d = h/oj = r (l + 8 2 H 2 /h 2 ) l/2 (see Fig. lA16> . In the limit that h 2 /B 2 H 2 <s 1, this can be solved to give 

h (ajBHr ) 1/2 , (A7) 
which can be combined with equation iA4l and the definition of H to give the following ODE for h: 

dh 3hro 2coBro 



dr rl+zl h 



= 0. (A8) 



Flows from Smooth Ionization Fronts 



19 




This can be solved directly by means of the substitution y — h and the integrating factor (rf, + z%) to give 



1 + 



1/2 



which may in turn be substituted back into equation (IA4> to give 



r 2\ 

A'(r„) = A'„|l + -J 



From equation (IA6> . the streamline angles are then given by 



\l/2 

, ■ . w 1 r ° 
tan ; = — . 



(A9) 



(A10) 



(All) 



Hence, the vertical scale height on the axis is given by ho = {a>/3/2) 1 z*, at which point the streamlines are vertical. We can 
now fix the value of ll>/3 by comparing this with our numerical results from the model simulations. Figure0shows that for zero 
curvature, h^/zh — 0.4, which implies a>/3 0.3. Moving away fr om the axis, the streamlines start to diverge slightly, reaching 
an asympotic angle of i x = tan _1 (w/2/3)^ 2 as r — > oo. In Figure |A"T71 this behavior of the streamline angle is compared with 
that seen in our numerical simulations of Model B, where it is found that co//3 0.6 gives a reasonable fit for r ( > < z,, which, 
combined with the estimate of co/3, yields a> 0.4 and /? 0.7. 

The validity of the analytical resu lts for ro greater than a few times z* is not to be expected, since h no longer satisfies the 
conditions stated after equation (IA3> above. Nonetheless, it is satisfying that the general characteristics of the photoevaporation 
flows seen in our numerical simulations can be reproduced from these simple physical arguments. 

It is interesting to consider how the calculation would be affected by the presence of dust grains in the ionized flow, which 
absorb part of the ionizing radiation. To first order, this may be treated by considering t he co lum n dens ity of gas between the 
ionizing star and the ionization front, which is proportional to Nh/ cos 8. From equations iA9\ and (IA10> it can be seen that this 
column density is a constant for all values of j*o. As a result, dust absorption has no differential effect on the structure of the flow 
in this approximation but merely reduces the effective ionizing luminosity of the star. 

B. ARE THE IONIZATION FRONTS D-CRITICAL? 

The traditional classification of ionization fronts into weak/critical/strong and R-type/D-type (iKahnl 1 9541 iGoldsworthvl 1 96 ll) 
was carried out for plane-parallel fronts. Although the fronts can be curved in our simulations, the ionization front thickness is 
always much less than the radius of curvature and so a plane-parallel approximation is justified for the front itself. It is therefore 
an interesting question whether or not our fronts are indeed D-critical, as is often assumed for photoevaporation flows. 

Two approaches are possible to the empirical classification of the ionization fronts in the simulations. One possibility is to look 
at the velocity of propagation of the front with respect to the upstream gas. The other is to look at the internal structure of the 
front, w hich is the method we investigate first. 

Fi gure lBT8l graphs various quantities as a function of ionization fraction, x. Figure lBT8h shows that almost all the gas heating 
in the front occu rs at v ery low values of x, so the ionization transition itself can be considered approximately isothermal. Fur- 
thermore, Figure lBT8t > shows that the particle flux is roughly constant through a great part of the front. At very low ionization 
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Fig. B 18. — Structure of the ionization front: various physical quantities as a function of ionization fraction, x. (a) Gas temperature, (b) Particle flux, (c) Ionized 
density, (d) Isothermal Mach number, with inset showing the same quantity on a logarithmic scale for small values of x. The different models are shown by: A 
(solid line), B (dashed line), and C (dotted line). The analytic curves of ionized density and Mach number for a steady, plane-parallel, isothermal, D-critical front 
are shown by thick gray lines. All velocities are given relative to the heating front. 



fractions, the flow is not steady because the shock propagates at a slightly faster speed than the ionization front. At high ioniza- 
tion fractions, the streamlines start to diverge and the flow is no longer plane parallel, resulting in a reduction of the particle flux. 
Nonetheless, the flow can be seen to be approximately steady and plane-paralle l over virtually the en tire front. 

In this case, one can apply the simple model developed in Appendix A of Hennev et al. (2005J). For a time-steady, plane- 
parallel, isothermal, D-critical front, where the critical point occurs at x — 1, then the results of that paper imply 

2l/2_Q _ U/2 

M ^ = n M/2 ( B1 > 

(1 + xyi 2 

and 

«i « rpr. (B2) 

1 - [(1 -x)/2] 1 ' 2 

These curves are shown in Figure lBToT : and d (thick gray lines), together with the same quantities for the three numerical models 
at an age of 78, 000 yr. The normalization of the theoretical n\ curve is arbitrary. The three models are very similar and also agree 
quite well with the analytical curve, except for near x — and x = 1 . In particular, the D-critical curve implies that the Mach 
number immediately after the heating front at x - should be M n = V2 - 1 0.414, which is close to that observed in the 
numerical models. This implies that the models are all close to D-critical, although they all pass through the sonic point before 
reaching x = 1, with the most divergent model (C) doing so at x ^ 0.9. Model A, which is the least divergent, shows a somewhat 
smaller initial jump in the Mach number, which may be evidence that it is weaker than D-critical. 

The numerical models all show a maximum in the ionized density between x = 0.6 and x = 0.8, similar to that predicted by the 
D-critical curve. At higher values of x, on the other hand, the ionized density decreases much faster than in the analytical model 
because of the divergence in the numerical models. Note, however, that the particle flux is constant in the vicinity of the maxima 
in rii, so divergence is not important there. 

The other test of the D-critical hypothesis is to look at the Mach number on the neutral side, that is, the Mach number of the 
gas flow relative to the heating front as it enters the front from the shocked neutral shell, M n . For a steady, plane-parallel front, 
there exists a maximum value, M n max = c n /2ci, which corresponds to a D-critical front. M n is about 0.06 for the models shown 
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Fig. B19. — Evolution with time of relative Mach numbers of the gas flow either side of the heating front, (a) Relative Mach number with which neutral gas 
enters the heating front, M n (black lines). Also shown (gray lines) is the critical Mach number for a steady flow, c n /2ci. (b) Relative Mach number with which 
gas leaves the heating front, Mh, calculated at x = 0.1. Gray line shows the analytical approximation for a D-critical front. 



in Fig. lB18fc f. as can be seen from the inset. The evolution with time of M n for the models is shown in Figure lB"T9b . which also 
shows (gray line) the theoretical maximum value, M a max , taking c; to be the sound speed at the sonic point in the flow. It can be 
seen that M n actually exceeds M n max everywhere except for the late-time evolution of Model A! This should not be possible for a 
steady flow but is probably due to the fact that the sound speed in the neutral shell, c n , varies spatially within the shell, according 
to the past history of the propagation speed of the shock front, U s , thus violating the steady-state assumption even though the 
flow from the shell through the ionization front is approximately steady. Hence, M a is not a suitable diagnostic for classifying 

the ionization front. 

The evolution with time of Mh is shown in Figure lB"T9t >, calculated at the point where the gas leaves the heati ng fro nt at x — 0. 1 . 
At this point in the flow, the gas has already heated up to 10 4 K, although it is still largel y ne utral, see Figure IbToI The analytic 
value for a D-critical front with a sonic point at x = 1 is shown by the gray line (see eq. dBU above). The numerical results for 
Model B can be seen to trace the analytical line almost exactly, indicating that the front is indeed D-critical in this model. For 
Model C, Mh slightly exceed s that analytical value, but this can be explained by the fact that the sonic point occurs nearer to 
x = 0.9 in this model (see Fig. lB18fc 0. For Model A, Mh falls increasingly below the D-critical value as the evolution progresses, 
indicating that the ionization front becomes increasingly weak-D. This implies that there should be no sonic point within the 
ionization front itself. This is consistent with the fact that the sonic point occurs much further from the heating front than in the 
other two models (see Fig.|5j:), in a region where the gas is fully ionized. 



